# -*- coding: utf-8 -*-
import math, random
import matplotlib.pyplot as plt

NMILE = 1852
W = 4*NMILE
H = 2*NMILE
D0 = 110
alpha = 1.5*math.pi/180
theta_beam = 120*math.pi/180

OV_MIN, OV_MAX = 0.10, 0.20

# ---------------------------
# 基础函数
# ---------------------------
def depth(x):
    return D0 - (x - W/2) * math.tan(alpha)

def coverage_width(x):
    return 2*depth(x)*math.tan(theta_beam/2)

def rotate_point(x,y,theta):
    """坐标旋转 (顺时针theta)"""
    u = x*math.cos(theta)+y*math.sin(theta)
    v = -x*math.sin(theta)+y*math.cos(theta)
    return u,v

def rotate_back(u,v,theta):
    """逆旋转回原坐标"""
    x = u*math.cos(theta)-v*math.sin(theta)
    y = u*math.sin(theta)+v*math.cos(theta)
    return x,y

# ---------------------------
# 目标函数
# ---------------------------
def calc(solution):
    theta, offsets = solution
    offsets = sorted(offsets)

    if theta<0: theta+=math.pi
    if theta>=math.pi: theta-=math.pi

    rect = [(0,0),(W,0),(W,H),(0,H)]
    uv = [rotate_point(x,y,theta) for (x,y) in rect]
    umin,umax = min(u for u,v in uv), max(u for u,v in uv)

    intervals=[]
    for u0 in offsets:
        L=coverage_width(W/2)
        intervals.append((u0-L/2,u0+L/2))
    if not intervals: return 1e12

    intervals.sort()
    merged=[]
    a,b=intervals[0]
    for c,d in intervals[1:]:
        if c<=b:
            b=max(b,d)
        else:
            merged.append((a,b))
            a,b=c,d
    merged.append((a,b))
    covered=sum(b-a for a,b in merged)
    uncovered=max(0,(umax-umin)-covered)

    penalty=0
    for i in range(1,len(offsets)):
        L=coverage_width(W/2)
        d=offsets[i]-offsets[i-1]
        ov=1-d/L
        if ov<OV_MIN: penalty+=(OV_MIN-ov)*1e7
        elif ov>OV_MAX: penalty+=(ov-OV_MAX)*1e6

    total_len=len(offsets)*H
    cost=total_len+uncovered*1e7+penalty
    return cost

# ---------------------------
# 邻域扰动
# ---------------------------
def perm(solution):
    theta,offsets=solution
    theta+=random.uniform(-2,2)*math.pi/180
    offsets=list(offsets)
    if random.random()<0.5 and len(offsets)>1:
        i=random.randrange(len(offsets))
        offsets[i]+=random.uniform(-50,50)
    else:
        if random.random()<0.5 or len(offsets)<2:
            offsets.append(random.uniform(-500,W+500))
        else:
            offsets.pop(random.randrange(len(offsets)))
    return (theta,offsets)

# ---------------------------
# 绘制方案
# ---------------------------
def plot_solution(solution, overlaps):
    theta,offsets=solution
    rect=[(0,0),(W,0),(W,H),(0,H)]
    xs,ys=zip(*rect+[rect[0]])
    plt.plot(xs,ys,'k-')

    for i,u0 in enumerate(sorted(offsets)):
        vmin,vmax=-H,H*2
        x1,y1=rotate_back(u0,vmin,theta)
        x2,y2=rotate_back(u0,vmax,theta)
        plt.plot([x1,x2],[y1,y2],'b-')
    plt.rcParams['font.sans-serif'] = ['SimHei']
    plt.axis('equal')
    plt.title(f"θ={math.degrees(theta):.1f}°, 线数={len(offsets)}")
    plt.show()

    # 重叠率分布图
    plt.figure(figsize=(8,4))
    plt.plot(overlaps,'o-')
    plt.axhline(10,color='g',linestyle='--',label='下限10%')
    plt.axhline(20,color='r',linestyle='--',label='上限20%')
    plt.xlabel("相邻测线对编号")
    plt.ylabel("重叠率 (%)")
    plt.title("相邻测线重叠率分布")
    plt.legend()
    plt.grid(True,alpha=0.3)
    plt.show()

# ---------------------------
# 重叠率分布计算
# ---------------------------
def compute_overlaps(solution):
    theta,offsets=solution
    offsets=sorted(offsets)
    L=coverage_width(W/2)
    overlaps=[]
    for i in range(1,len(offsets)):
        d=offsets[i]-offsets[i-1]
        ov=1-d/L
        overlaps.append(max(0,ov*100))
    return overlaps

# ---------------------------
# 模拟退火
# ---------------------------
def simulated_annealing():
    theta0=0.0
    spacing=coverage_width(W/2)*(1-0.15)
    offsets=[]
    u=-200
    while u<W+200:
        offsets.append(u)
        u+=spacing
    cur=(theta0,offsets)
    cur_cost=calc(cur)
    best, best_cost=cur,cur_cost
    T=1e5; alpha=0.999
    rec=[]
    for it in range(20000):
        cand=perm(cur)
        cost=calc(cand)
        if cost<cur_cost or random.random()<math.exp(-(cost-cur_cost)/max(T,1e-9)):
            cur,cur_cost=cand,cost
            if cost<best_cost:
                best,best_cost=cand,cost
                print(f"iter {it}: new best={best_cost:.1f}, θ={math.degrees(best[0]):.2f}, n={len(best[1])}")
        rec.append(cur_cost)
        T*=alpha
    return best,best_cost,rec

# ---------------------------
if __name__=="__main__":
    best,best_cost,rec=simulated_annealing()
    print("最优角度:",math.degrees(best[0]))
    print("测线数量:",len(best[1]))
    print("总代价:",best_cost)

    overlaps=compute_overlaps(best)
    print("相邻测线重叠率 (%):",["{:.1f}".format(o) for o in overlaps])

    plot_solution(best, overlaps)

    plt.plot(rec)
    plt.title("SA 收敛曲线")
    plt.xlabel("迭代")
    plt.ylabel("代价")
    plt.show()
